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Abstract. Iterative phase retrieval algorithms typically employ projections onto 
constraint subspaces to recover the unknown phases in the Fourier transform of an 
image, or, in the case of x-ray crystallography, the electron density of a molecule. For 
a general class of algorithms, where the basic iteration is specified by the difference 
map, solutions are associated with fixed points of the map, the attractive character 
of which determines the effectiveness of the algorithm. The behavior of the difference 
map near fixed points is controlled by the relative orientation of the tangent spaces of 
the two constraint subspaces employed by the map. Since the dimensionalities involved 
are always large in practical applications, it is appropriate to use random matrix theory 
ideas to analyze the average-case convergence at fixed points. Optimal values of the 
7 parameters of the difference map are found which differ somewhat from the values 
previously obtained on the assumption of orthogonal tangent spaces. 
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1. Introduction 

In two dimensional interferometric imaging, as well as diffraction from three dimensional 
crystals and non-periodic objects, the Fourier modulus of an object is sampled on a 
discrete and finite grid. To recover the object, additional a priori constraints (object 
support, positivity, atomicity, etc.) are imposed in order to "retrieve" the corresponding 
phases, without which the Fourier synthesis of the object is impossible. A natural setting 
for computations is the N dimensional complex Euclidean space E , where N is the 
number of pixels (or voxels) in the finite Fourier transform grid. The most widely 
used phase retrieval algorithms are formulated as maps E N — > E , which, when used 
iteratively, impose the Fourier modulus data as well as the a priori constraints in the 
course of recovering the object. 

The work reported here concerns the convergence of the phase retrieval algorithm 
based on the difference map[[L]: 

p i ► D(p) = p + (3(U 1 of 2 -U 2 o /i) (p) . (1) 

The action of D is to add to the current object, p G E N , a difference of projections, 111 
and n 2 . In general, the action of a projection on an object p is the minimal modification 
of p that restores a particular constraint, where minimality is specified by the standard 
Euclidean norm in E N . The two specific constraints considered below are the Fourier 
modulus constraint, implemented by lip, and a "fixed" or "atomic" support constraint, 
implemented by lis. The nonzero real parameter f3 serves as the step size of the 
iterations. Convergence of the difference map is crucially dependent on the maps f\ 
and f'2 with which the basic projections in (1) are composed: 

f l {p) = {l+ ll )U i {p)- ll p (2 = 1,2). (2) 

The structure of the maps fi is geometrical, taking p to a general point parametrized 
by 7i on the line determined by p and IT(p). A previous analysis of the convergence 
of the difference map, greatly simplified by the assumption of orthogonality of the two 
constraint subspaces near the solution, found 

7i Pt = (3«) 

7 2 ° pt = l/P , (36) 

as optimal parameter values. The main result reported here is a modification of 
these values for the more realistic case of constraint subspaces that are not perfectly 
orthogonal, locally, in the vicinity of the solution. 

The presentation is organized as follows. Section 2 reviews the relationship between 
the fixed points of the difference map and the solution to the associated phase retrieval 
problem. Explicit expressions for the two chief projections and their linearization are 
derived in Section 3. Convergence of the difference map for a general pair of linearized 
projections is studied in Section 4. Traces encountered in Section 4 are averaged in a 
standard random matrix ensemble in Section 5, and compared with averages over object 
ensembles appropriate for phase retrieval. Results are summarized in Section 6. 
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2. The difference map 

A special case of the difference map, the hybrid input-output map^, first appeared 
in the context of two dimensional phase retrieval with a fixed support constraint. The 
generality of the construction for an arbitrary pair of project ions [III |3|, and the flexibility 
in choosing the parameters 7t0, was only noticed recently. With IT = lis (support), 
n 2 = n F (Fourier modulus), the hybrid input-output map is obtained for the parameter 
values 71 = -1, 72 = 1/(3 Q. 

A phase retrieval solution p so1 exists if and only if the difference map has a fixed 
point p*. This is an immediate consequence of the definition of the difference map. At 
a fixed point p* = D(p*), the difference of projections in (p]) vanishes, and we have 

n 1 o/ 2 (p*) = n 2 o/ 1 (p*) : =p s ° 1 . (4) 

From this follows 

n^p 801 ) = n 2 (p so1 ) = p so1 , (5) 

showing that the object p so1 satisfies both constraints. Conversely, given p so1 , the set of 
fixed points is given by 

(IT o / 2 )-V o1 ) n (n 2 o /O-V 01 ) , (6) 

a nonempty set since it trivially contains p so1 . 

To recover an object using the difference map one begins with an initial, arbitrary 
object p(0), and generates iterates p{n) = _D n (p(0)) until the norm of the difference 

\\ p{n + i) _ p ( n )|| = ||0 (Hl o f 2 - II 2 o f\) (p(n))\\ (7) 

becomes suitably small. Then, from the approximate fixed point p* ~ p(n), the solution 
p so1 is obtained using @. There is no guarantee, that given an arbitrary starting point 
p(0), the iterates will converge on a fixed point. In particular, it may happen that the 
linearization of the difference map in the neighborhood of a fixed point has unstable 
directions. Avoiding this, and instead making fixed points maximally stable in an 
average sense, is the primary focus of the present work. 

The uniqueness of solutions is linked to the tightness of the constraints. In many 
applications the two constraint subspaces C\ and C 2 are unions of countably many 
submanifolds of E N and their stable intersection properties {i.e. solutions) are functions 
of their dimensionalities. Thus if dim(Ci) + dim(C 2 ) > N, then the set of solutions, 
C\ n C 2 , has positive dimension and there is no uniqueness. For phase retrieval to 
be well posed, the number of a priori constraints must be sufficiently large, and the 
dimensionality of the associated constraint subspace correspondingly small, so that only 
the empty intersection is stable. In this overconstrained situation a solution will reduce 
to a single point (or low dimensional submanifold in the case of symmetries) since a 
solution is known to exist at the outset. 
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3. Projections 



When the Fourier transform of an object is sampled on a regular grid (fundamentally 
unavoidable in crystallography, an artifact of CCD array design otherwise), the object 
can be treated as living on a periodic real-space grid. The fineness of the grid in real- 
space is determined by the range of the Fourier-space sampling. A general point of the 
periodic real-space grid has the form 

r = (m, r 2 , . . .) r i eZ m ./m i , (8) 

where the integers rrii give the ranges of the Fourier-space sampling for the different 
components of the Fourier-space basis, and Z m denotes the set of integers modulo m. 
Points in Fourier-space (relative to a basis which may be non-orthogonal) have the form 

q = (q h q 2 , . . .) q i e2-K r L mi . (9) 

In crystallography, r is called a "fractional coordinate" , while q is, up to a factor of 2tt, 
a triplet of Miller indices. The periodicity of the Fourier-space grid @ is an artifact of 
the regular, grid-like sampling in real-space; this artificial periodicity can be tolerated 
since the modulus of the Fourier transform of real objects becomes negligible at the 
extremes of the range. With points in real- and Fourier-space defined as above, the 
components of an arbitrary object in real-space, p r , and Fourier-space, p q , are related 
by the unitary transformation (discrete Fourier transform) with matrix elements 

(F)qr= exp it? -r , (10) 
where N = m\ m 2 • • • is the total number of real- or Fourier-space samples. 



3.1. Support projection 

Object support constraints can take two forms, depending on the application. In phase 
retrieval with non-periodic objects one normally has a priori knowledge of a fixed object 
support S, with the property p r = if r S. For example, relatively tight bounds 
on the set S can be derived from the object's autocorrelation as given directly by 
the Fourier transform of the squared Fourier modulus. In crystallography, on the other 
hand, it is the form of the support that supplies the necessary constraint. Given the 
number of atoms M, a valid support comprises a union of M very compact sets, say 
3x3x3 arrays of voxels. Although the locations of the individual atomic supports 
is unknown, the small size of the combined support of all M atoms is a very strong 
constraint. 

In the case of a fixed object support the corresponding projection operator is linear 
and has the following simple form: 

(Ils)rr/ = /]S ra 6 ar ' ■ (H) 

ses 

It is easily verified that given an arbitrary object p, the object II s (p) is the nearest object 
(in the Euclidean sense) having support S. The corresponding constraint subspace is a 
linear subspace of E N with dimensionality \S\. 
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The constraint subspace for an atomic support is, in contrast, a union of very many 
linear subspaces, each one corresponding to a different arrangement of atoms. Thus one 
is really working with a collection of supports S, where the linear projection operator 
(PI) applies to each S G S. The actual projection operator is nonlinear, returning the 
element of {IIs(p):S G S} which minimizes ||Hs(p) — p||. This complication, however, 
is not relevant to the local analysis we perform below: only one element of S contains 
the solution (given a choice of object origin and enantiomer) and the linear projection 
(PI) can be used for just that one. 



3.2. Fourier modulus projection 

We consider the situation where a large number of the object's Fourier moduli have 
been measured to high precision. In this case, the moduli F q in 

p s 9 ol = F q expi0 q , (12) 

for a large subset Qdata of Fourier samples q, can be treated as known quantities. The 
corresponding projection operator is nonlinear but simple when expressed in Fourier- 
space: 

' P 

Fq— 3 - for q G Q da ta and p q ^ 
~ \Pq\ 

n F (p q ) = < Fq for q G Q da ta and p q = ( 13 ) 

p q otherwise. 

On actual computing machinery, and random initialization, the arbitrary choice of phase 
made in the second case above never comes up in practice. In crystallography, Qdata 
never includes the origin and frequently also omits "reflections" q close to the origin. 
Since in all applications there is a systematic decay of moduli at large spatial frequencies, 
a significant fraction of the samples near the corners of the Fourier grid (those farthest 
from the origin) will lack measured values. There is usually no harm in approximating 
these small moduli by zero, since it is the absolute error that matters in Fourier synthesis. 
Geometrically, Fourier modulus projection in Fourier-space corresponds to projections 
of the complex numbers p q onto circles with prescribed radii (F q ), for all q G Qdata- A 
more refined approach, one which takes into account measurement errors and bounds 
on unmeasured moduli, would involve projections, respectively, onto annuli and disks. 
When considered together with support projection (|TT|), it is necessary to use the real- 
space form of Fourier modulus projection: 

n F = J^- 1 on F o^. (14) 



3. 3. Linearized projections 

To study convergence, the difference map and the projections from which it is built are 
linearized about the solution p so1 . In general, the linearization tt of a projection II is 
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,(,):=li m n ^' + t "'-^'. (15) 



defined by 



Since object support projection (|TT| ) is already linear, 

TTsfa) = n S (7?) . (16) 

Linearized Fourier modulus projection (0) is diagonal in Fourier-space, 

{§ - \ exp (i20 q ) C for q G Qdata and F q 7^ 
for q G Qdata and F q = (17) 

1 otherwise, 

where C is the complex conjugation operator. The second case above, although never 
strictly encountered with realistic objects, expresses the fact that the rank of the 
projection effectively decreases by one whenever \p q \ ^> F q (projection onto the origin, 
rather than a circle). This case is included because small (unmeasured) moduli at 
large spatial frequencies are ubiquitous in practice. Assuming for simplicity, however, 
that Qdata includes all N Fourier samples, and that all the moduli F q are nonzero, the 
real-space form (]TJj) of linearized Fourier modulus projection is given by 

Mw' = - 7^ exp i [20 q - qr • (r + r')} C . (18) 

Q 

To facilitate the comparison with random matrix theory, we express 7Tp in somewhat 
more compact notation. In the case of real-valued objects, where the operator C can 
be dropped, we have 

7r F = ^(l + U T U), (19) 

where the unitary matrix U is defined by 

(U) qr =-±=ex V i( ( j )q -q-r). (20) 



Since for real objects 0_ 9 = —<p q , the reality of U T U follows from 

(uy qr = (u)^ qr . (21) 

When representing complex- valued objects in a 2iV-dimensional real vector space, a 
general operator X + YC, where X and Y are ordinary matrices, has the block structure, 

Re(X + Y) lm(Y-X) 



(22) 



lm(X + Y) Re(X-Y) 

and the projection takes the form 

1 + Re(U T U) lm(U T U) 
lm(U T U) 1 - Re(U T U) 

The linearized difference map is defined by. 



D(p so1 + en) - p so1 

din) := lim {P - U P — , (24 

e^O e 
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and is obtained by replacing the projections Hi and II 2 in D by their linearized 
counterparts, 7i"i and ir 2 - The fixed points of d form a linear space that can be directly 
constructed from the kernels of 7Ti and 7r 2 as follows. From (^) and p* = p so1 + rj* we 
have 

TTi o [(1 + 72 )tt 2 - 72 ] (77*) =0 (25a) 
7r 2 o[(l + 7l ) 7 r 1 - 7l ](^) = 0. (256) 

By forming suitable linear combinations of fl25a|) and ( |25fe|) , and making use of the 
idempotency of tti and 7r 2 , one can show 

tti(O = (26) 

7r a (ifi = . (27) 

Thus, 

ker (d — 1) = ker 7Ti PI ker 7r 2 := ker 7Tj_ , (28) 

where ir± is the projection onto the orthogonal complement of the linear space of fixed 
points. Since it is the projection of the iterates d n (r)) by tc±_ that should converge to 
zero, we consider in the next Section the operator 

d± := 7Tj_ o d o . (29) 
Applying DeMorgan's law (for Euclidean subspaces) to (^|), we obtain 

ker (1 - tti) U ker (1 - tt 2 ) = ker (1 - tt ± ) . (30) 
By using ([28]) and ([30]) one easily verifies a set of identities needed in the next Section: 

7Ti O 7Tj_ = 7T_L O 7Tj = 7Tj (z = 1, 2) . (31) 

Finally, since we confine ourselves to the overconstrained case where 

ker (1 - tti) n ker (1 - vr 2 ) = {0} , (32) 

it follows from (|30|) that 

dim ker (1 — 7Ti) + dim ker (1 — 7T 2 ) = dim ker (1 — 7Tj_) , (33) 

or equivalently, 

Tr 7Ti + Tr 7r 2 = Tr 7Tj_ . (34) 
4. Convergence of the difference map 

With linearization, convergence to the space of fixed points corresponds to the vanishing 
of the norm of the iterates {d±) n (r]), where the map d±_ was defined in (F2"9"p. This in turn 
is guaranteed by a suitable bound on an induced matrix norm of d±. With the standard 
Euclidean norm on the vector space, the natural choice for the induced matrix norm is 
the spectral norm: 

IM_l|U := max\d±(r])\ . (35) 

|7?| = 1 
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There are two reasons, however, why the spectral norm may not be the appropriate 
choice. First, it presents a serious challenge to calculate, even in the limit iV — > oo. 
Second, although being able to achieve j| || oo < 1 by suitable choice of parameters 
guarantees convergence, it may not be true that minimizing this bound also maximizes 
the average rate of convergence. An alternative is suggested by the interpretation of 
(H^jJoo) 2 as the standard oo-norm applied to the eigenvalues of the matrix d±' o d±. 
In particular, the 1-norm of the (non-negative) eigenvalues has the advantage that it is 
easily calculated: 



Id i 



\2 ._ 



— Trrf_L f o d A 



(36) 



The norm called the Frobenius norm, gives the root-mean-square change in the 

Euclidean length of uniformly sampled unit vectors when acted on by d±. Although a 
bound on does not guarantee convergence, its minimization makes more sense in 

that it reflects a property of all the eigenvalues rather than just the top of the spectrum. 
This feature, and its tractable form, leads us to adopt the Frobenius norm as the vehicle 
for optimizing the difference map; the subscript 1 will subsequently be dropped. 
A straightforward calculation, making repeated use of (|3~1~D , gives 

\\d x \\ 2 = t ± + 2p [ 7l t 2 - 72*1 + (72 - 7i)*i2] + 

(3 2 [ifh + 72 2 tx + (2 + 2(71 + 72) - (71 - 7 2 ) 2 ) ti2 - 2(1 + 7l )(l + 72)^1212] 

where 



(37) 



(38a) 
(386) 
(38 c) 
(38 d) 
(38e) 

Since ( |3~?| ) is quadratic in the 7, and non-negative, it has a unique minimum where 
the partial derivatives <9(||dj_|| 2 )/<97i vanish. The optimal ji are given by the following 
expressions: 

{ti — t 12 )(t 2 — t 1212 ) — (*12 — ^1212) (*1 — 2ti2 + *1212)/3 



h 


= (Tr7n)/A^ 




t 2 


= (Trvr 2 )/iV 




t± 


= (TT7T ± )/N = t 1 


+ t 2 


tl2 


= (Tr 7ii H2) / N 




^1212 


= (Tr 7Tl O 7T2 O 7Tl 


O 7T 2 )/N 



opt 

7i 



opt 
72 



[tl{t 2 -t 12 )+t 
^12) (*1 ~~ *1212) — (tl2 



12(2*1212 — l>2) - > I.2I.2J 



t 2 )~t\ 
tl21 2 )(t 2 — 2t\ 2 



h 2 i 2 )P 



^1212] P 



[t 2 (ti — ti 2 ) + ^12(2^1212 - ti 

evaluated for these parameter values is 

|| cf_i_ || 2 = t± 

+ {(h-h 2 )(t 2 -t 12 )(2t 1212 -t 1 

+ [(tl2-tl212)(2tlt 2 -3ti 2 (tl 



(39 a) 
(396) 



The Frobenius norm 

12 



t 2 ) + 2{t 2 -t 1 )(t 12 -t 1212 ) 2 f3 



At 



12 



'-1212 



(ti + t 2 - 2t X2 )] P 2 }/ 



\t\t 2 — t\ 2 (t\ + t 2 — 2^1212) — ^1212} • 



(40) 
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When the linearized constraint subspaces are orthogonal, as was assumed in reference , 
the product traces t 12 and i m2 vanish and one recovers the result 7° pt = — 
72 Pt = 1/(3. Using fl38cj) in (|40|), one obtains ||d±|| 2 = for this case. This simple 
result could have been anticipated since for the given values of the 7, and orthogonal 
subspaces (tti o tt 2 = 0) one has the simplification d = 1 — -k\ — -k 2 = 1 — 7Tj_, and hence 
d ± = 0. 



5. Average traces of projection products 

Because our convergence measure ||g?j_|| 2 is linear in the traces ti 2 and £1212, an average 
of these traces with respect to a particular class of objects will provide the optimum 
average-case convergence. The projections that appear in these traces can be expressed 
as 7r = T PO, where P is a fixed canonical projection and O is an orthogonal 
matrix. The actual matrices O that apply to the projections considered earlier comprise 
subensembles in the space of orthogonal matrices. For tt s (object support) O is 
a permutation matrix, while for 7Tf (Fourier modulus) O belongs to a continuous 
family of orthogonal matrices defined implicitly by ([19]) or (^). Traces of products 
of 7Ti = Oi 1 P\0\ and n 2 = 2 T P 2 2 only depend on the relative orthogonal matrix 
O = Oi0 2 T . In the absence of detailed knowledge of the appropriate subensemble for 
O, we might consider the average over the complete, group invariant ensemble: 

(tia) := ^(Tr P 1 T P 2 0) (41a) 

:= ^(Tr P 1 T P 2 OP 1 T P 2 0) . (416) 

There are highly sophisticated methods 0, || for evaluating these standard random 
matrix averages which rely solely on symmetry. Here we quote the limiting forms for 
large N; details of the calculation are given in the Appendix: 

lim (t 12 ) = ht 2 (42a) 

N^oo 

lim (t 1212 ) = t\t 2 + t x t\ - t\t\ . (426) 

In the remainder of this Section these formulas are compared with explicit averages over 
actual object ensembles (and 7Ti = 7Ts, tt 2 = 7Tp). If the random matrix theory results 
( |42q| , [42 b\) are reproduced, then the problem of choosing optimal parameters is greatly 
simplified since (^7[) will then only depend on the dimensionalities of the constraint 
subspaces (ti and t 2 ). 

We begin by calculating the traces of the two projections. From (]TT|) and (|i~6l) we 

have 

t s = (Trrc s )/N=\S\/N:=a , (43) 

where a is the fraction of pixels or voxels in the support S. 

For traces involving the linearized Fourier modulus projection ( |18D we need to 
understand the action of the conjugation operator C. In the case of real- valued objects 
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all traces may be evaluated in a real vector space of dimension N, and C can be replaced 
by the identity. For complex-valued objects traces are evaluated in a 2iV-dimensional 
real vector space where the representation of a general operator X + YC has the block 
form (p2|). Normalizing the matrix trace of ( |22| ) by an extra factor of | gives us the 
identity 

Tr (X + YC) := Tr (ReX) , (44) 

where the trace on the right is an ordinary matrix trace. Another useful identity 
involving C is 

(X + YC)(W + ZC) = (XW + YZ*) + (XZ + YW*)C , (45) 

where X, Y, W, and Z are ordinary matrices. 

Use of identity ([44]) on ([18]) gives us immediately 



t F = (Trvr F )/iV = i (46) 
in the case of complex-valued objects. For real-valued objects we have 

(Trvr F )/iV = i-^^exp(i20 (J )E r exp(-i2q-r)] . (47) 



The inner sum in ( [4TD vanishes unless the components of q satisfy g.j = (mod 7rm 
and there are at most 4 (two dimensions) or 8 (three dimensions) such points. In more 
concrete terms, the second term in ( f4"T|) subtracts the number of samples q which do 
not possess a continuous phase when the object is real. This correction is negligible in 
the N —>■ oo limit and we recover the complex-object result ([46]). 

We consider next the first product trace computed earlier using random matrix 
theory. In the complex-object case, identities (|4~4] ) and ( ^3[ ) give us immediately 



t S F = (Tr7r s o7r F )/iV= °- , (48) 

in agreement with the random matrix result ( [42 cj) , even without averaging. For real- 
valued objects we have 

(Tr tts o vr F )/iV = | 1 J>P (i20 q ) S_ 2q , (49) 



where 



S<7 ' N 

ses 

defines the Fourier transform of the support's characteristic function. 

The second term in ( |4"9"P is a correction to the naive random matrix theory result 
and in general does not vanish. To check whether it vanishes in an average-case sense, 
we consider two object ensembles distinguished by the nature of their support. If the 
support is a known fixed set S, say a rectangle or disk in two dimensions, then S_ 2(? is a 
constant in the ensemble average of (|49|) and we can focus on the average of exp (i20 q ). 
The minimal a priori knowledge of the object in real-space is that the individual pixel 
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(or voxel) values p r in the support are independently distributed according to a known 
distribution. For this object ensemble it is straightforward to calculate the distribution 
of the Fourier transform components p q . Excepting a fraction of samples q which 
vanishes as iV — > oo (near q = or associated with facets in the support), the complex 
numbers p q are isotropic Gaussian random variables with mean zero and corrections 
that decay in the limit iV —>■ oo. Consequently, the average of the phases exp (i2(/> q ) in 
([|9|) vanishes in the same limit, as does the correction (since |£ q | < a, independent of 
N). 

In the case of an atomic support, although the actual support S is not known, one 
does know that 5* is the union of a given number M of very compact sets associated 
with atoms. In the equal-atom ensemble there is an approximate relationship 

S 9 ~ ~^Pl ' ( 51 ) 



where A(q) is a smooth real-valued function (form factor) independent of N. The 
approximation in (|51|) arises from shifts (by fractions of voxels) of the atomic centers 
relative to their individual supports. To bound the sum in fl4"9|) we replace each term 
by its modulus and use (|5l|) : 



1 

2iV 



eX P ( i2 0q) S -2q 



1 V |A(2g)| 



q 

The root-mean-square of the modulus \p2q\ can be calculated in the standard atomic 
ensemble where M atom centers are uniformly and independently distributed; the result, 
|p2q|rms = 0(^/M/N), (2q ^ 0) shows that the correction is negligible for N — ► oo and 
fixed resolution (M/N = const). We conclude that explicit averaging in the standard 
equal-atom ensemble, as in the fixed support ensemble, reproduces the random matrix 
theory result 



(*sf> = \ ■ (53) 

Finally, we turn to the trace of the product of four projection operators. For 
complex-valued objects we obtain 

^sfsf = (Tr 7T S o tt f o 7T S o tt f )/N = j + — ^expi2(0 q - (j) ql ) Y? q ,_ q , (54) 
and for real objects: 

^SFSF = (Tr 7T S O 7T F O 7T S O Tl F )/N = J - — ^ eX P Wq) S -2q 

q 

i-^expi2(0 g + (7 OSV< J - ( 55 ) 



AN 
q,q' 

We recognize the single sum in (p5| ) as the sum in P5|), which was already argued to 
be negligible when averaged over the standard object ensembles with fixed or atomic 
support. However, the averages of the double sums for the fixed support ensemble, both 
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in ( |54| ) and (|55|), can deviate from the random matrix theory result ( |42 b\ ) . In fact, the 
latter is recovered in both cases by retaining only the diagonal terms in the double sums 
(q = q' in (|54|), q = —q' in fl55|)). Since the contributions of the near-diagonal terms 
is not diminished in the N — > oo limit, there will be non- negligible corrections that 
depend on details of the ensemble (support shape, distribution of object values, etc.). 
The most we can conclude for the fixed support ensemble without these details is, 

(W> = j + 0(a 2 ) , (56) 

a result that may be useful for small object support. 

The neglect of non-diagonal terms in the double sum may be easier to justify 
rigorously for the atomic support ensemble. Considering only the case of real atoms 
(|55|), we first use fl5"T|) to write 

^-q'-q ~ \P-q'-q\ eX P (-^<Pq'+q) ■ (57) 

Now, since \p- q '- q \r ms = 0{M/N) for q + q' ^ 0, the non-diagonal terms have moduli 
which scale as N~ 2 (at constant resolution). To argue that the 0(N 2 ) non-diagonal 
terms in the double sum have a negligible ensemble average, we consider the phase of a 
general term: 

exp i2(0 q + <p ql - <p q ' +q ) . (58) 

Atomic ensemble averages of triplet phase invariants such as ( |58"D are studied at length 
in the statistical theory of structure factors in crystallography^]; they are found to 
vanish in the limit of many atoms, or equivalently, iV — > oo when the resolution is fixed. 
There are three exceptions to this general result, when q = 0, q' = 0, or q' + q = 0. 
The latter corresponds to the diagonal contribution in the double sum of (|55|) which we 
will retain. The two other cases lead to single sums whose moduli can be shown to be 
negligible by arguments identical to those used for the sum in fl52|). In summary, this 
non-rigorous argument leads us to believe that for the atomic support ensemble 

2 

(*sfsf> = \ + °j , (59) 
in agreement with the random matrix theory result 



6. Optimal parameter values 

The results of the previous Section show that random matrix averaging, of the traces 
of products of projections, reproduces for the most part the averages over the actual 
object ensembles encountered in standard phase retrieval applications (fixed and atomic 
support). This greatly simplifies the expressions for the optimal 7 parameters of the 
difference map (|39o| , |39l] ), since all the traces are explicit functions of just the traces 
ts = o~ and £f = h 

4 + (2 + 3)a + 3a 2 
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6 _ 2 g - (3(2 - 3a + a 2 ) 

7f = (3(4- a + a 2 ) • (606) 

Following the discussion of uniqueness in Section 2, phase retrieval with a support 
constraint is well posed when t s + t F < 1, which requires that the support fraction 
satisfies < a < |. As noted in Section 5, with the fixed support ensemble the average 



of tsFSF O —) ma y differ from the random matrix average by terms proportional to 



a 2 , and the coefficients of a 2 in the above formulae will change accordingly. For the 
atomic support ensemble the formulae as written are more reliable. Moreover, in the 
principal atomic support application, crystallography, there is a large dimensionality 
asymmetry of the two constraint subspaces as expressed by the statement (7 < ^ F° r 
this important application one may therefore take the a — > limit: 

7s Pt =° ~ (61«) 

7 opt^o3__^ (gi6) 



The Frobenius norm (P0|) in this limit is given by 

||d ± || 2<J =°^(3 + 2/3 + 3/3 2 ). (62) 



The numerical constants in (|6T£| ) and (^2|) can be traced to the numerical value of 
t-F = ~, the only trace with a finite value in the a — > limit. Recalling the discussion in 
Section 3.3 and the interpretation of tpN as the dimensionality of the Fourier modulus 
constraint subspace, the exact value ~ is certainly an oversimplification. As a first 
improvement in the evaluation of Tr tt^ one should exclude contributions from samples 
q for which the corresponding modulus F q is negligible. In crystallography, for example, 
only the samples within an ellipsoid about the origin have non-negligible F q , effectively 
reducing the rank of the Fourier modulus projection by a fraction given approximately 
(in three dimensions) by n/6 ~ 0.52; thus tp ~ 0.26. Given the uncertainty in t F we 
give formulas for the optimal 7 parameters for general tp: 

7s° Pt CT =° ~ (63a) 
7 o P t ^0 l+t F (l-0) m 
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Appendix 

Averages over products of matrix elements of N x N orthogonal matrices O (N > 1), 
such as 

can always be expressed as tensor products of Kronecker-5's having the correct invariance 
properties. The latter follow from the invariance of the group integration measure with 
respect to the change of variables O — > LOR, where L and R are independent orthogonal 
matrices. For the example above, the statement of invariance becomes 

44 X Z R i R i = X H > ( A - 2 ) 

with general solution 

X;{ = a5 ij 5 pq , (A.3) 

where a is an undetermined constant. By contracting the indices i and j, whereupon 
the left hand side of ( |A.3|) corresponds to ((0 T 0) pq ) = 6 pq , consistency requires 

a = l/N. (A.4) 

Using this method for the product of four matrix elements, 

XUri = {0\0\O k r O l s ) , (A.5) 



we find 



Kir's = ai(^5 kl 5 pg 6 rs + 8 tk V l 6 pr 6 qs + 5 il P% s 8 qr ) + (A.6) 
a 2 (5 ij 5 kl 5 pr 5 qs + ^5 kl 5 ps 5 qr + 5 ik V l 5 pq 5 rs + 
6 ik 6% s 5 qr + 6 u V k 5 pq 6- rs + 5 a V k 5 pr 5 qs ) . 

The reduction to two undetermined constants made use of permutation symmetries. 
Contracting indices k and / in ( |A.5|) gives our previous result for ( |A.1|) multiplied by 
5 rs , which can then be compared with the result of contracting the same pair of indices 
in (|A.6|) . Consistency yields two equation with the solutions 

N + l 

ai ~ (N-l)N(N + 2) [ ' 

a2 = (N-l)N(N + 2) • (A ' 8) 
The random matrix averages of the traces of products of projections are obtained 
by contracting the tensor expressions ( |A.3| ) and (|A.6|) with the projections P\ and P2: 

(Tr P 1 T P 2 0) = a(TiP l )(TiP 2 ) (A.9) 

(Tr P 1 T P 2 OP 1 T P 2 0) = a 1 (TrP 1 )(TrP 2 )(l + TrPx + Tr P 2 ) + (A.10) 

a 2 (TrP 1 )(TrP 2 )(3 + TrP x + TrP 2 + (Tr P 1 )(TrP 2 )) . 

After taking the limit — > 00 in (|A.10|) we obtain the results quoted in Section |. 
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